Calculation of the conductance of a graphene sheet using the Chalker-Coddington 

network model 
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The Chalker-Coddington network model (introduced originally as a model for percolation in the 
quantum Hall effect) is known to map onto the two-dimensional Dirac equation. Here we show how 
the network model can be used to solve a scattering problem in a weakly doped graphene sheet 
connected to heavily doped electron reservoirs. We develop a numerical procedure to calculate the 
scattering matrix with the aide of the network model. For numerical purposes, the advantage of 
the network model over the honeycomb lattice is that it eliminates intervalley scattering from the 
outset. We avoid the need to include the heavily doped regions in the network model (which would 
be computationally expensive), by means of an analytical relation between the transfer matrix 
through the weakly doped region and the scattering matrix between the electron reservoirs. We 
test the network algorithm by calculating the conductance of an electrostatically defined quantum 
point contact and comparing with the tight-binding model of graphene. We further calculate the 
conductance of a graphene sheet in the presence of disorder in the regime where intervalley scattering 
is suppressed. We find an increase in conductance that is consistent with previous studies. Unlike 
the tight-binding model, the network model does not require smooth potentials in order to avoid 
intervalley scattering. 

PACS numbers: 73.50.Td, 73.23.-b, 73.23.Ad, 73.63.-b 
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I. INTRODUCTION 

The low-energy and long-wave-length properties of 
conduction electrons in a carbon monolayer (graphene) 
are described by the two-dimensional Dirac equation^. In 
one-dimensional geometries this partial differential equa- 
tion can be solved analytically, but fully two-dimensional 
problems typically require a discretization to permit a 
numerical solution. The tight-binding model on the hon- 
eycomb lattice of carbon atoms provides the most obvi- 
ous and physically motivated discretization^. The band 
structure of a honeycomb lattice has two valleys, cou- 
pled by potential variations on the scale of the lattice 
constant. Smooth potentials are needed if one seeks to 
avoid inter-valley scattering and obtain the properties of 
a single valley. 

Discrete representations of the Dirac equation that 
eliminate from the outset the coupling to a second valley 
may provide a more efficient way to isolate the single- 
valley properties. Alternative tight-binding models^^^ 
have been introduced for that purpose. One method 
of discretization which has received much attention is 
the network model, originally introduced by Chalker and 
Coddington as a model for percolation in the quantum 
Hall effect 7 . Ho and Chalker 8 showed how a solution of 
this model can be mapped onto an cigcnstate of the Dirac 
equation, and this mapping has proven to be an efficient 
way to study the localization of Dirac fermions 9 . 

The recently developed capability to do transport mea- 
surements in graphene^ has renewed the interest in the 
network modeUi and also raises some questions which 
have not been considered before. The specific issue that 
we address in this paper is how to introduce metallic 



contacts in the network model of graphene. Metallic con- 
tacts are introduced in the Dirac equation by means of 
a downward potential step of magnitude [Too- The limit 
Uqo — > oo is taken at the end of the calculation. (It 
is an essential difference with the Schrodinger equation 
that an infinite potential step produces a finite contact 
resistance in the Dirac equation.) This phenomenologi- 
cal model of metallic leads, introduced in Ref. [HI, is now 
commonly used because 1) it is analytically tractable, 
2) it introduces no free parameter, and 3) it agrees well 
with more microscopic model a 13 i 14 . A direct implemen- 
tation of such a metallic contact in the network model is 
problematic because the mapping onto the Dirac equa- 
tion breaks down in the limit Uqc — > oo. Here we show 
how this difficulty can be circumvented. 

To summarize then, there is a need to develop numer- 
ical methods for Dirac fermions in graphene when the 
potential landscape does not allow analytical solutions. 
If one implements a method based on the honeycomb 
lattice of graphene, intervalley scattering is present, un- 
less the potential is smooth on the scale of the lattice. 
Smooth potential landscapes are experimentally relevant, 
but computationally expensive, because they require dis- 
cretization with a large mesh. It is therefore preferable 
to develop a numerical method that eliminates intervalley 
scattering from the outset. The known correspondence 
between the Chalker-Coddington network model and the 
Dirac equation provides such a method, as we show in 
this paper. The key technical result of our work is an 
analytical method to include heavily doped reservoirs. 
(Including these reservoirs numerically would have been 
prohibitively expensive, computationally.) 

In Sees. |TT] and IIIII we summarize the basic equations 
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FIG. 1: Top panel: Schematic of a graphene sheet contacted 
by two electrodes. A voltage source drives a current through 
the sheet. The bottom panel shows the potential profile 
V(x,y) for fixed y. 



that we will need, first regarding the Dirac equation and 
then regarding the network model. Our key technical re- 
sult in Sec. IIVI is a relationship between the scattering 
problems for the Dirac equation in the limit Uoo oo 
and for the network model at Uoo = 0. We test the 
method in Sec. |V] by calculating the conductance of 
an electrostatically defined constriction (quantum point 
contact) in a graphene sheet. We also study the effect of 
disorder on conductance. We confirm the results of previ- 
ous studie a 15 i 16 ' 17 i 18 that smooth disorder (that does not 
cause intervalley scattering) enhances the conductivity of 
undoped graphene. We conclude in Sec. IVil 



II. FORMULATION OF THE SCATTERING 
PROBLEM 

A. Scattering Matrix 

A scattering formulation of electrical conduction 
through a graphene sheet was given in Ref. [T^- We 
summarize the basic equations. The geometry, shown in 
Fig. [T] consists of a weakly doped graphene sheet (length 
L and width W) connected to heavily doped graphene 
leads. A single valley has the Dirac Hamiltonian 

H = vcr ■ [p — eA(r)] + V(r) + cr z (j,(r), (2.1) 

where A(r) is the magnetic vector potential, V(r) is the 
electrostatic potential, and /i(r) is a substrate-induced 
mass term. The vector cr = [a x ,a y ) contains the stan- 
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We assume that the fields A, V, and /i are smooth on 
the scale of the lattice constant, so that the valleys are 
uncoupled. 

In the heavily doped leads (for x < and x > L) we set 
V(r) — —Uoo and take the limit Uoo — > oo. For simplicity 
we set [l — in the leads and we also assume that the 
magnetic field is zero in the leads (so A is constant there). 
The Dirac equation 

H^ = E^ (2.3) 

has to be solved subject to boundary conditions on the 
wave function ^(r) at y — and y — W . We will con- 
sider two types of boundary conditions which mix neither 
valleys nor transverse modes. The first is the periodic 
boundary condition \I/ | y=0 = w . The second is the 
infinite-mass boundary condition 24 

*l»=o = *l w =o . ^y=w = *\ v =w ■ ( 2 - 4 ) 

We consider a scattering state ty n that has unit inci- 
dent current from the left in mode n and zero incident 
current from the right. (The quantum number n labels 
transverse modes.) In the leads ^ n has the form 

*»(r) = xUv) e lknX + J2 r ™ e ~ ZkmX > x < °' 

m 

(2.5a) 

# n (r) = ]T t mn xtn(v) e^~ L \ X > L. (2.5b) 

III 

We have introduced transmission and reflection ampli- 
tudes t mn and r mn and the longitudinal component k n 
of the wave vector of mode n. The right-propagating 
component in mode n has a spinor \n an d the left- 
propagating component has a spinor %n ■ 

In the limit Uoo ^ oo, the form of the scattering 
state in the leads can be simplified considerably. The n- 
dependence of k n can be neglected, since k n ~ Uoa/fov — > 
oo as Uoo — * oo. The number Noo — UooW/hv of 
propagating modes in the leads can be taken infinitely 
large. When Noo — * oo, the choice of boundary condi- 
tion in the leads (not in the sample) becomes irrelevant 
and we choose periodic boundary conditions in the leads 
for simplicity. Modes that are responsible for transport 
through the weakly doped sample have transverse mo- 
menta \q n \ <C Uoo- The corresponding spinors \n are 

^=^^(±0- qn=2 ^ (2 - 6) 

with n = 0, ±1, ±2, . . .. While it is important not to 
neglect the finiteness of q n in the phase factor exp(iq n y) 
of these modes, the spinor structure is proportional to 
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(1,±1) independent of n, because q n /U a 
the orthogonality relation 



0. We note 



(2.7) 



We also note that the definition of Xniv) ensures that 
each scattering state ^ n carries unit incident current. 

In a similar way, we can define a scattering state in- 
cident from the right in mode n with transmission and 
reflection amplitudes t' mn and r' mn . The transmission and 
reflection amplitudes constitute the scattering matrix 



S = 



r t' 
t r' 



(2- 



which is a unitary matrix that determines transport prop- 
erties. For example, the conductance G follows from the 
Landauer formula 



4e 2 ■ 4e 2 
G= ^Trttt = ^Tr t't'\ 
h h 



(2.9) 



where the factor of 4 accounts for spin and valley degen- 
eracies. 



B. Transfer matrix 

The information contained in the scattering matrix S 
can equivalently be represented by the transfer matrix 
T. While the scattering matrix relates outgoing waves to 
incoming waves, the transfer matrix relates waves at the 
right, 

= E KXn(yy ak " (x ' L \ x > L, (2.10) 

n, a 

to waves at the left, 

*i(r)=£<^(y)e iCTfe " x , x<0. (2.11) 



where E z is a matrix in the space of modes with entries 
(S 2 ) m = 5 mn a z that are themselves 2x2 matrices. In 
terms of the transfer matrix the Landauer formula 
can be written as 



4e 2 
h 



(2.15) 



C. Real-space formulation 

In order to make contact with the network model, it is 
convenient to change from the basis of transverse modes 
(labeled by the quantum number n) to a real space basis 
(labeled by the transverse coordinate y). The real space 
transfer matrix X y _ y i is defined by 



IV 



dy' X V!y ,*(0,y'), (2.16) 



where &(x, y) is any solution of the Dirac equation 
at a given energy E. The kernel X y>y i is a 2 x 2 ma- 
trix, acting on the spinor \&. Because the integral (|2.16|) 
extends only over the weakly doped region, X does not 
depend on the potential in the leads. 

In view of the orthogonality relation (|2.7[) the real- 
space transfer matrix X is related to the transfer matrix 
T defined in the basis of modes in the leads by a projec- 
tion Onto Xmi 



m.n 



W rW 



dy / dy'x?n(yyX y , y ,Xn(y')- (2-17) 



We now substitute the explicit form of xZ from Eq. (|2.6[) . 
The integrals over y and y' in Eq. (|2.17|) amount to a 
Fourier transform, 

X m , n = — dy dy 1 e-^ y X y . y ,e^y' . (2.18) 
" Jo Jo 



The relation takes the form 



b a = T <T ' cr 'a cr ' 

m / j m,n ' 



(2.12) 



The four blocks T a,a of the transfer matrix are related 
to the transmission and reflection matrices by 



r = -(T— ) _1 T- 4 
t = T++ - T+- (T~ 

a = (T-y\ 

r' = T+- (T—)~\ 



T" 



(2.13a) 
(2.13b) 
(2.13c) 
(2.13d) 



Unitarity of S implies for T the current conservation re- 
lation 



T- 1 = E ? T f S, 



(2.14) 



From Eq. (|2.17[) we conclude that the 2x2 matrix 
structure of the transfer matrix, 



rp+ + 

rp / J m.n m,n 

^m,n I rp j- rp— — 

m.n m.n 



(2.19) 



is related to the 2x2 matrix structure of the real-space 
transfer matrix by a Hadamard transformation: 



(2.20) 



(The unitary and Hermitian matrix Tt is called the 
Hadamard matrix.) In view of Eq. (|2.14p . the current 
conservation relation for X reads 



X 1 — Ti x X'H x , (Sa,) 
where we used TLa z TL = a x . 



,„,„ = S m>n a x , (2.21) 
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FIG. 2: Square lattice (dots), with circulating current loops 
that form the network model. The loops are coupled to near- 
est neighbors at the black rectangles. The lattice vectors ai 
and a 2 (each of length \/2Z) are indicated. 



% (m+ l,n) 




(m, n — 1) 



FIG. 3: Segment of the network of Fig. [5] with the wave am- 
plitudes Zm]n and scattering matrices „ indicated. 



III. FORMULATION OF THE NETWORK 
MODEL 

The Chalker-Coddington network modefr^ was origi- 
nally introduced in order to analyze the localization tran- 
sition in the quantum Hall effect. Our interest in this 
model stems from the fact that it is known to map onto 
the two-dimensional Dirac equation^ We briefly recall 
how the network model is defined and how the mapping 
to the Dirac equation works. We consider the square lat- 
tice shown in Fig. [2] with lattice constant V21 and lattice 
vectors 



ai = l(x + y), a 2 = l(y-x) 



(3.1) 



The matrix elements of s+ n and s m n are arranged 
such that 




7 (i) 

-* m n — 1 

7(3) 



(2) 



e^m+l.n 



e*£U 
e**&. 



S m,n I „(3) 



(3.2a) 



7(2) 



S m,n I ^(4) 



m,n — 1 



(3.2b) 



The integers (m, n) label the lattice site r OT , n = tocii + 
na,2- With each site is associated a single current loop 
circling the site without enclosing any neighboring sites, 
say clockwise if viewed from the positive z axis. The 
radii of these loops are expanded until states associated 
with nearest neighboring sites overlap. At these points 
of overlap, states on adjacent loops can scatter into each 
other. 



As illustrated in Fig. [31 four current amplitudes Z m>n , 



m,n). These 



k = 1, . . . , 4 are associated with each site 
are amplitudes incident upon points of overlap, ordered 
clockwise, starting from the point of overlap with site 
(to+1, n). Each incident wave amplitude Zmjn has picked 
up a phase 4>m]n since the previous point of overlap. With 
the point of overlap between loop (to, n) and (to + 1, n) 
is associated a 2 x 2 scattering matrix s+ n , while s~ n is 
associated with the point of overlap between (to, n) and 
(to, n — 1). 



Ho and Chalker 8 showed how this model can be 
mapped onto the Dirac equation for two-dimensional 
fermions. Firstly, one parametrizes the scattering ma- 
trices n in terms of Pauli matrices <t^, 



Bin Q + P m ,n^j <7» + COS Q + /3 m ,„) (T x , 

(3.3a) 

s m,n = COS (J + /3 m , n J <7 a + sill Q + /3 m ,„) CT a . 

(3.3b) 



(The same matrix of coefficients /3„ l; „ is used for s+ „ 
and s~ .) For given fields V(r), A(r), and fi(r) in the 
Dirac equation, the mapping then dictates a correspond- 
ing choice of parameters in the network model, namely 
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and P m ,n have to satisfy 8 

^X>^„ - [B-K(r„ lB )] 

k=l 

6 (1) -6 (3) el 

6 {A) ~6 {2} el 
— = A v {r m , n ) — , 



l 

hv 



hv ' 



(3.4a) 

(3.4b) 

(3.4c) 
(3.4d) 



With this choice of parameters there is an approximate 
equality between a solution ^(r) of the Dirac equation 
and the current amplitudes of the network model, 



order in this quantity.) This would imply that very large 
networks are required for an accurate representation of 
the graphene strip. 

It turns out, however, that it is not necessary to model 
the heavily doped leads explicitly in the network model, 
as we now demonstrate. We define the real-space transfer 
matrix Y as the matrix that relates Z^ and Z^> at the 
right edge of the network to and Z^> at the left 
edge of the network. The left edge (x — 0) lies at sites 
(n, n) with n = 0, 1, 2, . . . , N — 1. The right edge at 
x = L = 2 Ml lies at sites (n + M,n- M). The real- 
space transfer matrix Y relates 



•(i) 

n+M,n-M 



(3) 

n+M,n-M 



N-l 



( 

v I n'n' 
n'=0 \ n',n' 



(4.1) 



*(r m ,„) » g 



7 (3) 



Q = 



= -( l 

V2 V 1 



(3.5) 



The accuracy of the approximation is improved by mak- 
ing the lattice constant \/2l smaller and smaller. 

As mentioned in Sec. [Til we will be considering two 
types of boundary conditions at y — and y — W in 
the sample region < x < L. The periodic boundary 
condition is realized in the network model by putting the 
square lattice on a cylinder of circumference W — 2NI 
oriented along the a;-axis. The infinite-mass boundary 
condition is realized 8 - by terminating the square lattice 
at y — and y — W and adjusting the scattering phases 
along the edge. The edge y = lies at sites (n, —n) and 
the edge y — W lies at sites (N — 1 + n, N — 1 — n). 
As shown in App. [A] for sites (n, — n) Eq. (|3.2|) must be 
replaced with 

7 (4) _ j,(3) 7 (3) _ 7 (2) ,„ fi x 

while for sites (N + n, N — n) it must be replaced with 

7 (2) _ 7 (1) 7 (4) _ y{l) 

N+n,N — n ~ ^N+n.jV-n' ^N+n.N-n ~ ^N+n^N-n' 

(3.7) 



IV. CORRESPONDENCE BETWEEN 
SCATTERING MATRICES OF DIRAC 
EQUATION AND NETWORK MODEL 

In this section we combine the known results summa- 
rized in the previous two sections to construct the scat- 
tering matrix S of a graphene strip with heavily doped 
leads from a solution of the network model. This con- 
struction does not immediately follow from the corre- 
spondence (|3.5[) because the limit Uoo — > oo of heavily 
doped leads still needs to be taken. At first glance it 
would seem that, in order to preserve the correspondence 
between the network model and the Dirac equation, wc 
must simultaneously take the limit / — > so that UoqI/Tiv 
remains small. (The correspondence between the net- 
work model and the Dirac equation is correct only to first 



We define the Fourier transform 

N-l N-l 



Y* m ,in ^EE e- 2il «™ m 'Y m ,, n ,e 2il ^ n ', (4.2) 

rci 1 — n'= 

with q n = 2irn/W. 

In view of the relation (|3 . 5|) between the Dirac wave 
function ^> and the network amplitudes Z^\ Z^\ the 
real space transfer matrix X of the Dirac equation is 
related to Y by a unitary transformation, 



(4.3) 



We can now use the relation (|2.20p between X and the 
transfer matrix T to obtain 



where we have used 



1 

i 



Y„. 



1 

-i 



ng = 



1 o 

i 



(4.4) 



(4.5) 



From Eq. (|4.4|) it follows that the lower right blocks of 
T and Y are equal: T~ 
Landauer formula (|2.15[) gives 



y qmt q n - Substitution into the 



4e 2 
G = ^-Tr 
h 



y—\ Y — 



(4.6) 



The Landauer formula applied to the network model 
thus gives the conductance of the corresponding graphene 
sheet connected to heavily doped leads. For later use, we 
note the current conservation relation for Y, which fol- 
lows from Eqs. (|2~14"1) and (|4~4]) 



(4.7) 



y- 1 = s z y+s 2 



NUMERICAL SOLUTION 



In this section we test the accuracy and efficiency of 
the solution of a scattering problem in graphene by means 
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W/L 



FIG. 4: Transmission probability of a clean graphene sheet, 
at energy E — 7. 85 Tiv/L as a function of transverse wave 
number q. The solid line is the result (|5.3[l from the Dirac 
equation, while the open circles were numerically calculated 
using the network model with periodic boundary conditions 
(when q = 2im/W). The discretization parameter of the 
network was e — El/Tiv = 0.28. 



FIG. 5: Conductivity a — G x L/W at E — ior a, clean 
graphene sheet as a function of the aspect ratio. The data 
points were calculated from the network model for fixed L = 
40 Z with periodic boundary conditions (circles) and infinite 
mass boundary conditions (squares) in the weakly doped re- 
gion. The solid lines are the result from the Dirac equation. 
The dashed line indicates the limiting value ah/Ae 2 = 1/n for 
short wide samples. 



of the network model. As explained in Sec. IIVI we need 
to calculate the real space transfer matrix Y through 
the weakly doped region. The conductance of the corre- 
sponding graphene sample then follows from Eq. (|4.6p . 

We calculate the real-space transfer matrix recursively 
by adding slices to the network and multiplying the trans- 
fer matrices of individual slices. Since a multiplication of 
transfer matrices is numerically unstable we stabilize the 
algorithm as explained in App. [Bj We limit the numer- 
ical investigation in this section to the case A(r) = 0, 
/i(r) = where only the electrostatic potential V(r) is 
non-zero. 

We have found that the efficiency of the algorithm 
can be improved by using the fact that, according to 
Eq. (|3.4p , there is some arbitrariness in the choice of the 
phases 4>^ x \ . . . , (j)^ . For A(r) = and /i(r) = 0, one 
choice of the phases could be 

^]n = [E-V{ma 1 +na 2 )]l/2, fc = l,...,4. (5.1) 

Another choice is 

= A = [E- V(r m , n )} I, 4> {2) = (4) = 0. 

The correspondence (|3.5[) between the network model 
and the Dirac equation holds for both choices of the 
phases, however the corrections for finite I are smaller 
for choice (|5.2|) . More precisely, as shown in App. [Cl if 
</>( 2 ) and <j)W are zero, the network model does not con- 
tain corrections to the Dirac equation of order d r Vl. 

Let us first consider the analytically solvable case of a 
clean graphene sheet that is obtained by setting V = 
in the weakly doped region. The Dirac equation gives 



transmission probabilitie: 



T(E,q) = 



,12 




(5.3a) 
(5.3b) 



For periodic boundary conditions the transverse wave 
vector is discretized as q n = 2nn/W, with n = 
0, ±1, ±2, ... 

In Fig. 0] we compare Eq. (|5.3[) to the results from the 
network model for periodic boundary conditions in the 
weakly doped region. The small parameter that controls 
the accuracy of the correspondence is e = El/hv. We 
find excellent agreement for a relatively large e ~ 0.3. 

Fig. [5] shows the conductivity 



W h 



(5.4) 



at the Dirac point [E = 0) as a function of the aspect 
ratio W/L. We do the calculation both for periodic and 
infinite mass boundary conditions in the weakly doped 
region. (In the latter case q n — (n + ^)ir/W with n = 
0, 1, 2, . . .) Again we see excellent agreement with the 
analytical results from the Dirac equatioiJ^. 

We now apply the network model to a case that can- 
not be solved analytically, because it involves inter-mode 
scattering. We take the electrostatic potential landscape 
shown in Fig. [SI which produces a narrow constriction 
or quantum point contact of width D and length L c . In 
the weakly doped region, of length L, electrons have an 
energy Ep measured from the Dirac point. The barrier 
potential is tuned so that electron transport through the 
barrier takes place at the Dirac point, where all waves are 
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V = 



FIG. 6: Potential landscape V(x, y) that produces a quantum 
point contact. The Fermi energy Ef is indicated. 
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FIG. 7: Conductance through the constriction of Fig. [6] 
as a function of the width of the opening in the constric- 
tion. The solid line was obtained using the network model, 
while the dashed line was obtained using the tight-binding 
model of graphene. We used parameters W = 35 Hv/Ef, 
L c = 8.7 hv/ Ef ■ For the network model we set the length 
of the weakly doped region to L — 49 Tiv/Ef and used a lat- 
tice constant \/2l = 0.24 Hv/Ef, while in the tight-binding 
calculation we used a lattice constant 0.17 Hv/Ef- 



evanescent. As the constriction is widened, the number 
of modes at a given energy that propagates through the 
opening increases. For fixed Ep, this should lead to steps 
in the conductance as a function of opening width, at in- 
tervals of roughly n/Ep. The steps are smooth because 
the current can also tunnel through the barrier. 

We have calculated the conductance with the network 
model (solid curve in Fig. [7]) and using the tight-binding 
model of graphene (dashed curve). In the tight-binding 
calculation we did not connect heavily doped leads to the 
weakly doped region. This does not affect the results, as 
long as L ^> L c . 

Both calculations show a smooth sequence of steps in 
the conductance. The agreement is reasonably good, but 
not as good as in the previous cases. This can be under- 
stood since the tight-binding model of graphene is only 
equivalent to the Dirac equation on long length-scales. 

The final numerical study that we report on in this pa- 
per involves transport at the Dirac point through a disor- 
dered potential landscape. Recent experimental studies^ 



have observed electron and hole puddles in undoped 
graphene. The correlation length of the potential is larger 
than the lattice constant, hence intervalley scattering is 
weak. We are therefore in the regime of applicability of 
the network model (which eliminates intervalley scatter- 
ing from the outset). 




L 



FIG. 8: Illustration of the model of electron and hole pud- 
dles in a graphene strip that we have studied. The sample 
is divided into tiles. The value of the potential on a tile is a 
constant, here indicated in gray-scale, uniformly distributed 
between — Vmax and Vmax. The potential on different tiles 
is uncorrelated. We choose a mesh for the network such that 
each tile has size 10 I x 10 1, where the network lattice constant 
is y/2l. 



To model the electron and hole puddles, we devide the 
sample into an array of square tiles (Fig |SJ), where each 
tile has size 10 Z x 10 Z, y/2l being the lattice constant of 
the network model. The electrostatic potential is con- 
stant on a single tile, but uncorrelated with the potential 
on the other tiles. We take the values of the potential 
on any given tile to be a random variable uniformly dis- 
tributed between — V max and Vnax- To make contact with 
previous studie a 15 > 16 , we quantify the disorder strength 
by the dimensionless number 



K 



1 



{hv) 2 



dr' (V(r)V(r')) • 



(5.5) 



(The average (V(r)) is zero.) With tiles of dimen- 
sion 10 1 x 10 1, the relation between Kq and Vmax is 
Kq = 1 00(V max l/hv) 2 /3 and the network model faithfully 
represents the Dirac equation for values up to Ko ~ 10. 
We use a sample with aspect ratio W/L = 5 and average 
over 100 disorder realizations. We repeat the calculation 
for two different sample sizes namely W = 5L = 3001 
and W — 5L — 450L The calculation is performed for 
transport at energy E = 0, i.e. the Dirac point of a 
clean, undoped sample. In Fig. [5] we show the aver- 
age conductance. Remakably enough the conductance 
increases with increasing disorder strength. This is con- 
sistent with the results obtained in Refs. HEHIOIlI- 
The effect should not depend on the shape of the tiles in 
our model for the disorder. We have therefore repeated 
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FIG. 9: Conductivity a — G L/W averaged over 100 disorder 
realizations versus disorder strength Ko at the Dirac point 
E = 0. The circles are for samples of size 60 1 x 300 1 while 
squares are for samples of size 90 / x 450 1. The statistical 
error is of the order of the size of the data points. The dotted 
line indicates the ballistic limit G L/W = 4e 2 /nh. 



a honeycomb lattice. In agreement with the existing 
literatur o 15 i 16 ' 17 i 18 we have found that disorder that is 
smooth on the scale of the graphene lattice constant en- 
hances conductivity at the Dirac point. The absence of 
intervalley scattering in the network model may prove 
useful for the study of these and other single- valley prop- 
erties. 
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APPENDIX A: INFINITE-MASS BOUNDARY 
CONDITION FOR THE NETWORK MODEL 



the calculation with rhombic instead of square tiles. We 
find deviations of less than 5%. 

The increase in conductance is explained by the non- 
zero density of states at the Dirac point that is induced by 
the disorder, together with the absence of back-scattering 
for Dirac electrons. While we do not make a detailed 
study of the dependence of conductance on sample size 
(at fixed aspect ratio), we note that the conductance 
of larger samples (squares in Fig. [5]) is larger than the 
conductance of the smaller samples (circles in Fig. [S]). 
This is consistent with the scaling behavior found in 
Refs. QOEl 



VI. CONCLUSION 

In conclusion, we have shown how the Chalker- 
Coddington network model can be used to solve a scatter- 
ing problem in a weakly doped graphene sheet between 
heavily doped electron reservoirs (which model the metal- 
lic contacts) . The method is particularly useful when the 
scattering problem does not allow an analytical solution, 
so that a numerical solution is required. The network 
model eliminates intervalley scattering from the outset. 
Thus, with a given mesh size, a larger graphene sam- 
ple can be modeled with the network model than with 
methods based on the honeycomb lattice. The key tech- 
nical result of our work is that an infinitely high potential 
step at the contacts can be implemented analytically by 
a unitary transformation of the real-space transfer ma- 
trix, without having to adjust the lattice constant of the 
network model to the small values needed to accommo- 
date the small wave length in the contacts. We have 
demonstrated that the algorithm provides an accuracy 
and efficiency comparable to the tight-binding model on 



In this appendix we consider the boundary condition 
imposed on the Dirac equation by termination of the net- 
work along a straight edge. We consider the eight orienta- 
tions shown in Fig. ITD1 which have the shortest periodicity 
along the edge. Since we want to discuss the long wave- 
length limit, each edge needs to be much longer than 
the lattice constant y/2l. (In this respect the figure with 
its relatively short edges is only schematic.) The orienta- 
tions are defined by the vector n(a) — —x sin a + y cos a, 
a = j7r/4, j = 1 , . . . , 8 which is perpendicular to the edge 
and points outwards. 

We wish to impose the infinite mass boundary 
condition^ 

*ed ge = [n(a) x z ] ■ a * e dge 

= (a x cos a + a y sin a) * e dge ( Al ) 

on the Dirac wavefunction at the edge. In view of the 
correspondence (|3.5p between the Dirac equation and the 
network model, Eq. (|A1 1) implies the boundary condition 

I Z ( 3 ) J = {-Vx sma + a z cos a) I ^ (3) I 

V / edge V / edge 

(A2) 

on the network amplitudes. 

Away from the edge, the network amplitudes obey the 
equations ()3.2j) . For [i, A, V, and E all equal to zero 
(Dirac point) these reduce to 

( 7 {2) \ ( 7 {1) \ 

7 (4) = n \ 7 m . ( A3a ) 

( Sr 1 ) =n ( J> ) • (A3b) 

We can eliminate the amplitudes Z^> and to arrive 
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FIG. 10: Network of circulating current loops, as in Fig. [2] 
but now terminated with straight edges. The letters a, b, . . . 
label the orientation of the edge. 



at the equations 
7(1) = - 



ri+1 



,(1) 



7(3) , 7 (3) 



Z (3) 



There are 

7 (3) >> 



771 — 1,71 — 1 

, 7(3) 



m+l,ra+lj 



7 (3) 



^ m+l,n ' m,n — lj 

two linearly independent 

.(i) 7 (3) 



(A4a) 

(A4b) 
solutions 

m,n,Z^ n ) (X (1,0) and (Z£jn,Zm,n) 

When the network is truncated along an edge, the bulk 
equations (IA4j) do not hold for the amplitudes along 
the edge. We seek the modified equations that impose 
the boundary condition (|A2[) up to corrections of order 
{E-V)l/hv. 

The edge orientation a was previously considered by 
Ho and Chalker— . We consider here all four independent 
orientations a, b, c, and d. The other four orientations 
a', c', and d! are obtained by a symmetry relation. 

Edge a is constructed by removing all sites (to, n) with 
n> m. (See Fig. 1111 ) This means that the network am- 

(3) 

plitudes Zm rn are prevented from scattering into the non- 



existent amplitudes Z, 



(2) 



belonging to the removed 



sites (to — 1 , m). Similarly, the amplitudes Zm,m are pre- 
vented from scattering into the non-existent amplitudes 



,(3) 

trices s 



To do this one must modify the scattering ma- 



m— l,m 



so that Z^r^jYi 



can only scatter into Z, 



(4) 



I I 

I I 

I • I 

I (to, to + 1) 

I 



(m, TO ~r i- J 

f Z (2) Z (l) 

< — *L 



z^ z< 4 ' 



m,m+l 



(to, to) 



1 



Z (3) 



°m— l,m 1 

1 • 1 

' (m — 1, to) I 

I J z( 2 



FIG. 11: Network amplitudes at an edge with orientation a. 
The dashed current loops are removed. 



is replaced by 



(A6) 



The solution (zl'^zi^) cx (1,-1) indeed satisfies the 
infinite mass boundary condition (|A2|) with a = tt/2. 



(to,1) I 

J Z<$ Z^ 1 

A *L- ► 



7* 

Z& Z^iVS~ x 



(m,0) 



,(4) 



and s m m+l so that Z m , m 
a consequence, for n — to + 1 Eq. 



We eliminate Z^ 2 ^ and i?( 4 ) to arrive at Eq. (|A4 
n < to and Eq. (lA4b[) for n — m. Eq. (|A4a[) for n 



FIG. 12: Edge with orientation 6. 



Edge b is constructed by removing all sites (m, n) with 
n > 0. (See Fig. [l^]) This means that the network 
amplitudes Z^ are prevented from scattering into the 

(3) 



can only scatter into Z$ m . As non-existent amplitudes Z^ belongingtothe removed 



3J) is replaced by 

Zm]m- (A5) 



sites (m,l). For n = 1, we replace Eq. (|A3b[) by 



(4) 
m,0" 



(A7) 



for 

= TO 



If we now eliminate the amplitudes Z^ 2 ) and Z^ we find 
that Eq. Kil is still valid for all n < 0. For n = 0, 
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Eq. (|A4b|) still holds, while Eq. (|A4a[) is changed to 



yi 1 ) _ | y. 1 1 _ y'-" 



(1) 



.(3) 



(A8) 



The solution (zi'l, z£.l) T oc (1,1 — \/2) satisfies the 
infinite mass boundary condition (|A2[) with a = 7r/4. 



Z (4)| 
I 

(m, —to + V) j | ij« 

s 4 



• I 

(m, — to + 1) 



I 
I 

I • I 

I (m + 1, — m) | 

J 



Z< 3 > 



m, — m+ 



1 Z« 



(m, —to) 



Z (3) 
Z (2) 



FIG. 13: Edge with orientation c. 



Next, we consider edge c, which results from the re- 
moval of all sites (m, n) with to > — n. (See Fig. OJ) In 

this case, s m _ m+1 must be modified to prevent Z^_ m 

from scattering into Z^_ m+X . Furthermore, _ m 

must be modified to prevent Z^ _ m from scattering into 

Z { Si,-m- For n 



-m + 1 we replace Eq. (IA3|) by 



,(2) 

J m,—rn 



,(i) 



J ra.—ra 



,(4) 

J m,—m ' 



(A9) 



and eliminate i^ 2 ) and Z( 4 ) to verify that the boundary 
condition holds. 

The condition (|A9|) modifies three of the equations 
El: 



K3) 



ra,— ra — 1 



_ 7 (3) 
Tn— 1,— m ra.—ra. I ) 



l 

2 V ra,— 7ii— l m— 1,— m— 1 



-Z 



(3) 



ra. — ra} ' 

(i) 



(AlOa) 



(AlOb) 



ra — 1, — m— 1 



V2Z 



(1) N 
ra. — ra) 



(AlOc) 



For to < — rt — 1 Eq. (|A4[) holds without modification and 
Eq. (|A4bj) also holds for to = — n — 1 . The solution 

Zi]n<-m = ^ff- m = constant, Z<?>, =0 (All ) 

implies (Zm}n, Z^ n ) oc (1,0) for to < — n, which satisfies 
the infinite mass boundary condition (|A2|) with a = 0. 

Edge d results from the removal of all sites (n, to) with 
to > 0. (See Fig. |T4j) We must modify Sq m such that 




FIG. 14: Edge with orientation d 



%om does not scatter into Z^ m . To do this we replace 



Eq. (|A3a[) for sites (0, to) by 



,(2) 



(1) 



0,m 0,77i ' 

We again eliminate Z< 2 ) and 2J( 4 ) to arrive at 



(A12) 



7 (i) 



1 ^/2> (1) 7 (3 M 



(3) 
0,m 



(Al3a) 



) _ Z (D 



Z 



(3) 



m — l,m— 1 ~ 0,m.—l / ' 

(A13b) 



while for to < Eq. (|A4[) still holds. The solution 



(Z 



(1) 7 (3) 
m,n, ^m,n 



OC 



(1 , \/2 — 1 ) obeys the infinite mass bound- 



ary condition (|A2p with a = — 7r/4, as required. 

This completes the boundary conditions for the four 
orientations a, 6, c, and <i. The orientations a', 6', c', and 
d' are obtained by the following symmetry: The network 
model is left invariant by a n rotation in coordinate space 
(which takes r to — r) together with the application of 
(7y in spmor space (which takes to -iZ& and Z^ 
toiZ«). 



APPENDIX B: STABLE METHOD OF 
MULTIPLICATION OF TRANSFER MATRICES 

To construct the transfer matrix of a conductor one 
can divide it into slices, compute the transfer matrix of 
each slice, and multiply the individual transfer matrices. 
This recursive construction is numerically unstable, be- 
cause products of transfer matrices contain exponentially 
growing eigenvalues which overwhelm the small eigen- 
values relevant for transport properties. Chalker and 
Coddingtoni used an orthogonalisation metho d 20 ' 21 to 
calculate the small eigenvalues in a numerically stable 
way. To obtain both eigenvalues and eigenfunctions we 
employ an alternative metho d 16 ' 22 : Using the condition 
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of current conservation, the product of transfer matrices 
can be converted into a composition of unitary matrices, 
involving only eigenvalues of unit absolute value. 

We briefly outline how the method works for the real 
space transfer matrices Y of the network model, defined 
by Eq. (14.1 1) . For the recursive construction it is conve- 
nient to rewrite this definition as 



Z 



(i) 

m+L,m—L 



N-l 



7(3) 



"£y(l,l>) 7 



n=0 



7« 

J n-\-L' ,n—L' 
7 (3) 

J n-\-L',7i — L' 



(Bl) 

The numbers L, L' are integers, so that Y(L, V) is the 
transfer matrix from x' = 2L'l to x — 2LI. The compo- 
sition law for transfer matrices is matrix multiplication, 



Y(L,0) = Y(L,L-1)Y(L -1,0), 



(B2) 



with initial condition Y(0, 0) = identity matrix. 

The unstable matrix multiplication may be stabilized 
with the help of the condition Y~ x = Y, Z Y^Y> Z of current 
conservation (see Sec. lIVI) . Because of this condition, the 
matrix U constructed from Y by 



Y = 



a b 
c d 



O U = 



-dr x c 
- bd- 1 



d- 1 

c bd- 1 



(B3) 



is a unitary matrix (U 1 = W). Matrix multiplication 
of Y's induces a nonlinear composition of Z7's, 



defined by 



Oi h 
ci di 



Y 1 Y 2 ^U 1 <g)U 2 , 



a 2 h\ _ ( as h 
c 2 d 2 ) \c 3 d 3 



(B4) 



(B5) 



The algorithm now works as follows: Multiply a num- 
ber of transfer matrices and stop well before numerical 
overflow would occur. Transform this transfer matrix 
into a unitary matrix according to Eq. (|B3|) . Continue 
with the next sequence of transfer matrices, convert to a 
unitary matrix and convolute with the previous unitary 
matrix. At the end, we may transform back from U to 
Y by the inverse of relation (|B3|) 



U = 



A B 
C D 



& Y = 



C-DB- 1 A DB- 1 
-B~ X A B- 1 



(B7) 



In practice this final transformation is unnecessary. Ac- 
cording to Eq. (|B3|) the upper- right block of U is d" 1 = 
(Y -- ) -1 , which is all we need to calculate the conduc- 
tance using the Landauer formula (|4.6p . 
APPENDIX C: OPTIMAL CHOICE OF PHASES 
IN THE NETWORK MODEL 



In Sec.|V]we noted that the same long-wavelength cor- 
respondence between the Dirac equation and the network 
model can be obtained for different choices of the phases 
(j>m,n- Among these choices, the choice (15.2 [) avoids cor- 
rections of order d r V I to the Dirac equation. Here we 
show why. 



For n = A = Eq. ([3~4]) reduces to 

4>£]n = 4>t]n = (1 - «K 
( t > rn,n = 4 > m,n = a ^m.m 



(CI a) 
(Clb) 
(Clc) 



a 3 = a\ + &i(l — a2(ii) _1 a2Ci, (B6a) where we have defined the dimensionless quantity e m ,n = 

b b (1 ad )~ 1 b (B6b) [E — V{r m>n )\l/?iv. The parameter a can be chosen arbi- 

^ &2 1 ' 2 ' ^ ' trarily. We wish to show that the choice a — is optimal. 

c 3 = c 2(l - d ia2 y ci, (B6c) We substitute Eq. (|5^a|) into Eq. (pHb)) of Sec.QTIl with 

^3 = d 2 + £2(1 — dia 2 )~ 1 dib 2 . (B6d) this parametrization, and obtain 



zp - 
z (3 \ 



e ie m 



ni. a 



7 (3) 

" / m+l,n 



-ia(en 



-i-e m .„WD 



( 7^> - 7^ 

\^m-l,n-l ^ro.ra-l 



(C2a) 

(cab) 



Now we expand in e m>n , keeping terms to first order, and 
take Z^ and Z^ to be functions defined for all r and 
smooth on the scale of the lattice. From Eq. (|C2 [) we 



then obtain 



= [E + a z p x + <T x p y - V(r) 



a/V(r + o 2 ) - V(r) V(r + a 2 ) - V(r) 
2 I V(r) - V(r - a 2 ) V(r - a 2 ) - V{r) 



) 

(C3) 
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After transforming to * = Q(Z^\ Z (3) ) T , with Q as in 
Eq. (|3.5[) . the first term on the r.h.s. of Eq. (|C3|) becomes 
the desired Dirac equation. If we choose a^O then the 



potential V has to be smooth on the scale of the lattice, 
for the second term to be negligible in comparison with 
the first. We conclude that a = is the optimal choice. 
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